Cargamos los arhivos contenidos en “infoClin_all_dataset.RData” y “CRC_stageII.RData” con los contenidos clínicos y de CRC.
# Cargamos el archivo infoClin_all_dataset.RData 16/03/2024
dataClin <- load("data/infoClin_all_dataset.RData")
# Cargamos el archivo CRC_stageII.RData 19/03/2024 dataCRC <-
# load('data/CRC_stageII.RData')
# Cargamos el archivo exp_mat.RData 02/04/2024
dataCRC <- load("data/exp_mat.RData")
GSE Files
Platform GPL570: [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array
Fusionamos los datos en un único dataframe al que le añadimos la columna grupo con los datos stage y msi.
Factorizamos las variables categóricas y ordenamos por ID tanto los datos clínicos como los de la matriz de expresión.
Creamos un DataFrame para asignar cada ID a su estudio clínico. A cada estudio le asignamos un color para poder identificarlo fácilmente en los estudios gráficos.
# 11/04/2023 Utilizamos el paquete 'dplyr'
library(dplyr)
# Función para preparar el data frame con el estudio por ID. También
# asignaremos un color a cada estudio
prepare_df <- function(df, study_name, color) {
df %>%
select(ID) %>%
mutate(study = study_name, color = color)
}
# Lista de estudios con sus colores asociados
studies_info <- list(clx = "red", GSE39582 = "blue", GSE14333 = "green", GSE13294 = "orange",
GSE17536 = "yellow", TCGA = "violet")
# Creamos una lista de de los estudios utilizando la función prepare_df
df_studies_list <- lapply(names(studies_info), function(x) {
prepare_df(get(paste0(x, "_infoClin")), x, studies_info[[x]])
})
# Juntamos los df's de la lista se studies en uno solo y lo ordenamos por 'ID'
df_studies <- bind_rows(df_studies_list) %>%
arrange(ID)
# Fusionamos y ordenamos todos los datos clínicos generando 'df_dataClin'
df_dataClin <- bind_rows(lapply(names(studies_info), function(x) get(paste0(x, "_infoClin")))) %>%
arrange(ID)
# Cramos una nueva columna 'grupo' a los 'df_dataClin' con los datos de msi y
# stage
df_dataClin$grupo <- with(df_dataClin, paste(stage, msi_imputed, sep = "_"))
# Convertimos a factores
df_dataClin <- df_dataClin %>%
mutate(across(c(msi_imputed, cms, stage, grupo), factor))
# Añadimos las columnas 'study' y 'color' a df_dataClin
df_dataClin <- left_join(df_dataClin, select(df_studies, ID, study, color), by = "ID")
Pasamos la matriz de expresión a data_frame y la ordenamos por “ID”
# Pasamos la matriz de expresión a data frame 26/03/2024
df_exmat <- as.data.frame(ex)
# Ordenamos la matriz de expresión por nombre de columnas 26/03/2024
df_exmat <- df_exmat[, order(names(df_exmat))]
Resumen exploratorio de los datos:
# Hacemos un Summary con el paquete skimr 18/03/2024 install.packages('skimr')
require(skimr)
skim(df_dataClin)
| Name | df_dataClin |
| Number of rows | 1079 |
| Number of columns | 7 |
| _______________________ | |
| Column type frequency: | |
| character | 3 |
| factor | 4 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| ID | 0 | 1 | 7 | 12 | 0 | 1079 | 0 |
| study | 0 | 1 | 3 | 8 | 0 | 6 | 0 |
| color | 0 | 1 | 3 | 6 | 0 | 6 | 0 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| msi_imputed | 0 | 1.00 | FALSE | 2 | MSS: 840, MSI: 239 |
| cms | 158 | 0.85 | FALSE | 4 | CMS: 329, CMS: 262, CMS: 194, CMS: 136 |
| stage | 0 | 1.00 | FALSE | 2 | II: 684, III: 395 |
| grupo | 0 | 1.00 | FALSE | 4 | MSS: 511, MSS: 329, MSI: 173, MSI: 66 |
Estadísticas y cálculos para la generación de los gráficos de datos clínicos:
# Estadísticas desglosadas 18/03/2024 Frecuencias
msi_frecuencias <- table(df_dataClin$msi_imputed)
cms_frecuencias <- table(df_dataClin$cms)
stage_frecuencias <- table(df_dataClin$stage)
# Valores NA
na_msi <- sum(is.na(df_dataClin$msi_imputed))
na_cms <- sum(is.na(df_dataClin$cms))
na_stage <- sum(is.na(df_dataClin$stage))
# Contamos cuantos MSI y MSS tenemos por cada estadio, II y III 19/03/2024
sII_MSI <- sum(df_dataClin$msi_imputed == "MSI" & df_dataClin$stage == "II")
sII_MSS <- sum(df_dataClin$msi_imputed == "MSS" & df_dataClin$stage == "II")
sIII_MSI <- sum(df_dataClin$msi_imputed == "MSI" & df_dataClin$stage == "III")
sIII_MSS <- sum(df_dataClin$msi_imputed == "MSS" & df_dataClin$stage == "III")
# Crear un dataframe con estos valores
stage_MSI_MSS_df <- data.frame(Stage_MSI_MSS = c("II_MSI", "II_MSS", "III_MSI", "III_MSS"),
Freq = c(sII_MSI, sII_MSS, sIII_MSI, sIII_MSS))
# Cargamos la librería "ggplot"
require(ggplot2)
# Pasamos a DataFrame msi
df_msi_frecuencias <- as.data.frame(msi_frecuencias)
df_msi_frecuencias
Var1 Freq
1 MSI 239
2 MSS 840
# Añadimos los porcentajes
# Calculamos el total de frecuencias
total_frecuencias <- sum(df_msi_frecuencias$Freq)
# Añadimos una columna de porcentajes al dataframe
df_msi_frecuencias$Porcentaje <- (df_msi_frecuencias$Freq / total_frecuencias) * 100
# Diagrama de Barras de MSI/MSS
p1 <- ggplot(df_msi_frecuencias, aes(x = Var1, y = Freq)) +
geom_bar(stat = "identity", fill = c("salmon", "seashell3")) +
geom_text(aes(label = Freq), vjust = -0.5, size = 3.5) + # Frecuencias arriba de cada barra
# Porcentajes en la mitad + o -
geom_text(aes(label = paste0(round(Porcentaje, 1), "%")), vjust = 5, size = 3.5) +
labs(title = "MSI/MSS", x = "", y = "Frecuencia") +
theme_minimal()
p1
# Pasamos a DataFrame cms
df_cms_frecuencias <- as.data.frame(cms_frecuencias)
df_cms_frecuencias
Var1 Freq
1 CMS1 194
2 CMS2 329
3 CMS3 136
4 CMS4 262
# Añadimos los porcentajes
# Calculamos el total de frecuencias
total_frecuencias <- sum(df_cms_frecuencias$Freq)
# Añadimos una columna de porcentajes al dataframe
df_cms_frecuencias$Porcentaje <- (df_cms_frecuencias$Freq / total_frecuencias) * 100
# Diagrama de Barras de CMS
p2 <- ggplot(df_cms_frecuencias, aes(x = Var1, y = Freq, fill = Var1)) +
geom_bar(stat = "identity") +
geom_text(aes(label = Freq), vjust = -0.3, size = 3.5) + # Frecuencias arriba de cada barra
# Porcentajes en la mitad + o -
geom_text(aes(label = paste0(round(Porcentaje, 1), "%")), vjust = 5, size = 3.5) +
scale_fill_manual(values = c("salmon", "seashell3", "turquoise3", "slategray")) +
labs(title = "CMS", x = "", y = "Frecuencia") +
theme_minimal() +
theme(legend.position = "none") # Omitir la leyenda
p2
# Pasamos a DataFrame stage
df_stage_frecuencias <- as.data.frame(stage_frecuencias)
df_stage_frecuencias
Var1 Freq
1 II 684
2 III 395
# Añadimos los porcentajes
# Calculamos el total de frecuencias
total_frecuencias <- sum(df_stage_frecuencias$Freq)
# Añadimos una columna de porcentajes al dataframe
df_stage_frecuencias$Porcentaje <- (df_stage_frecuencias$Freq / total_frecuencias) * 100
# Diagrama de Barras de ESTADIO
p3 <- ggplot(df_stage_frecuencias, aes(x = Var1, y = Freq)) +
geom_bar(stat = "identity", fill = c("salmon", "seashell3")) +
geom_text(aes(label = Freq), vjust = -0.5, size = 3.5) + # Frecuencias arriba de cada barra
# Porcentajes en la mitad + o -
geom_text(aes(label = paste0(round(Porcentaje, 1), "%")), vjust = 5, size = 3.5) +
labs(title = "STAGE", x = "", y = "Frecuencia") +
theme_minimal()
p3
# Stage vs MSI/MSS DataFrame
stage_MSI_MSS_df
Stage_MSI_MSS Freq
1 II_MSI 173
2 II_MSS 511
3 III_MSI 66
4 III_MSS 329
# Añadimos los porcentajes
# Calculamos el total de frecuencias
total_frecuencias <- sum(stage_MSI_MSS_df$Freq)
# Añadimos una columna de porcentajes al dataframe
stage_MSI_MSS_df$Porcentaje <- (stage_MSI_MSS_df$Freq / total_frecuencias) * 100
# Diagrama de Barras de Stage vs MSI/MSS
p4 <- ggplot(stage_MSI_MSS_df, aes(x = Stage_MSI_MSS, y = Freq)) +
geom_bar(stat = "identity", fill = c("salmon", "seashell3",
"turquoise3", "slategray")) +
geom_text(aes(label = Freq), vjust = -0.5, size = 3.5) + # Frecuencias arriba de cada barra
# Porcentajes en la mitad + o -
geom_text(aes(label = paste0(round(Porcentaje, 1), "%")), vjust = 2, size = 3.5) +
labs(title = "STAGE vs MSI/MSS", x = "", y = "Frecuencia") +
theme_minimal()
p4
# Múltiple Layout con ggplot
require(patchwork)
# En un grid 2x2
p1 + p2 + p3 + p4 + plot_layout(ncol = 2, nrow = 2)
Para el Análisis de Componentes Principales (PCA), normalizamos y
escalamos (estandarizamos) la matriz de expresión. En nuestro caso, los
datos ya están normalizamos (en principio con \(log_2\)) por lo que tan solo los
escalamos.
La estandarización se realiza restando la media de cada variable y
dividiendo el resultado por la desviación estándar de esa variable.
\[
Z = \frac{x-\mu}{\sigma}
\] Utilizamos la función scale para escalar. A
continuación realizamos el PCA con la transpuesta de la matriz de
expresión dado que queremos un análisis de las muestras en función de
los genes y no al revés.
IMP: al hacer el escalado antes de la transposición
estamos estandarizando por genes y no por muestras.
# Escalamos/normalizamos la matriz de expresión. También centra. 28/04/2024 Al
# hacer el escalado antes de la transposición estamos escalando por genes y no
# por muestras.
df_exmat_scaled <- scale(df_exmat)
# Hacemos el PCA. Hay que hacer la transpuesta de la matriz de expresión (las
# columnas representan genes y las filas representen muestras), así los
# componentes principales reflejan la variabilidad entre muestras. 28/04/2024
# Na haría falta center = TRUE ya que al escalar ya centralizamos.
pca_result <- prcomp(t(df_exmat_scaled), center = TRUE, scale. = FALSE)
RESUMEN PCA
# Obtenemos los valores propios (la varianza explicada por cada componente)
varianza <- pca_result$sdev^2
# Calculamos la proporción de varianza explicada por cada componente
prop_varianza <- varianza/sum(varianza)
# Calculamos la proporción acumulativa de varianza explicada
prop_acumulada <- cumsum(prop_varianza)
# Crear un dataframe con esta información para las primeras 10 PC
resumen_pca <- data.frame(SD = pca_result$sdev[1:10], Varianza = prop_varianza[1:10],
VarAcumulada = prop_acumulada[1:10])
print(resumen_pca)
SD Varianza VarAcumulada
1 13.721545 0.08197744 0.08197744
2 11.124090 0.05387873 0.13585617
3 10.811433 0.05089263 0.18674880
4 9.600035 0.04012675 0.22687555
5 6.928531 0.02090117 0.24777672
6 6.605919 0.01900005 0.26677677
7 6.062676 0.01600358 0.28278035
8 5.642600 0.01386267 0.29664302
9 5.370890 0.01255975 0.30920276
10 5.178365 0.01167545 0.32087821
La baja varianza explicada por los primeros dos componentes, PC1
(8,2%) y PC2 (5.4%), sugiere que los datos tienen una distribución de
varianza muy dispersa y que no hay unas pocas direcciones principales de
variabilidad que dominen.
Esto es coherente con el hecho de que estamos tratando con datos
biológicos, en concreto de CRC, que suelen ser de alta dimensionalidad
con múltiples factores contribuyendo a la variabilidad general. Es
decir, no hay patrones dominantes de variabilidad que puedan ser
fácilmente capturados por los primeros componentes principales.
CARGAS DE GENES EN PCA
Para entender mejor qué genes están contribuyendo a la separación observada, podemos analizar las cargas de los componentes principales para ver qué genes tienen mayores pesos en los componentes que diferencian los grupos.
# Extracción de cargas de componentes principales
cargas_pca <- pca_result$rotation
# Mostramos las cargas de los primeros dos componentes
primer_componente <- cargas_pca[, 1]
segundo_componente <- cargas_pca[, 2]
PC1
# Extraemos los top 10 genes y sus cargas de PC1
top_genes_pc1 <- head(sort(primer_componente, decreasing = TRUE), 10)
# Convertimos a data frame para visualización de PC1
top_genes_pc1_df <- data.frame(Gen = names(top_genes_pc1), Carga = top_genes_pc1)
# Ordenamos el dataframe por carga de manera descendente
top_genes_pc1_df <- top_genes_pc1_df[order(top_genes_pc1_df$Carga, decreasing = TRUE),
]
# Ajustamos los factores para que el orden de los genes sea según las cargas
top_genes_pc1_df$Gen <- factor(top_genes_pc1_df$Gen, levels = top_genes_pc1_df$Gen)
# Creamos el gráfico de barras
plot_pc1 <- ggplot(top_genes_pc1_df, aes(x = Gen, y = Carga, fill = Gen)) + geom_bar(stat = "identity",
position = "dodge") + theme_minimal() + theme(axis.text.x = element_text(angle = 90,
hjust = 1, vjust = 0.5)) + labs(title = "Top 10 Cargas de Genes en PC1", x = "Gen",
y = "Carga") + scale_color_gradient(low = "peachpuff", high = "peachpuff4")
# Mostramos los top 10 genes de PC1
print(top_genes_pc1)
SFRP2 GAS1 COL10A1 CCL18 CYP1B1 SPP1 GPNMB
0.06248246 0.04743630 0.04494118 0.04386073 0.04165178 0.03902618 0.03896303
GREM1 CCDC80 POSTN
0.03860923 0.03848662 0.03844061
PC2
# Extraemos los top 10 genes y sus cargas de PC2
top_genes_pc2 <- head(sort(segundo_componente, decreasing = TRUE), 10)
# Convertimos a data frame para visualización de PC2
top_genes_pc2_df <- data.frame(Gen = names(top_genes_pc2), Carga = top_genes_pc2)
# Ordenamos el dataframe por carga de manera descendente
top_genes_pc2_df <- top_genes_pc2_df[order(top_genes_pc2_df$Carga, decreasing = TRUE),
]
# Ajustamos los factores para que el orden de los genes sea según las cargas
top_genes_pc2_df$Gen <- factor(top_genes_pc2_df$Gen, levels = top_genes_pc2_df$Gen)
# Creamos el gráfico de barras
plot_pc2 <- ggplot(top_genes_pc2_df, aes(x = Gen, y = Carga, fill = Gen)) + geom_bar(stat = "identity",
position = "dodge") + theme_minimal() + theme(axis.text.x = element_text(angle = 90,
hjust = 1, vjust = 0.5)) + labs(title = "Top 10 Cargas de Genes en PC2", x = "Gen",
y = "Carga") + scale_color_gradient(low = "peachpuff", high = "peachpuff4")
# Mostramos los top 10 genes de PC1
print(top_genes_pc2)
REG4 ZIC2 REG1A TRIM7 CTSE TCN1 AGR3
0.07056229 0.05049736 0.04890627 0.04815804 0.04715623 0.04230700 0.04170686
RAB27B SDR16C5 PLA2G2A
0.04125286 0.04028539 0.03911398
# Mostramos gráficamente las cargas del top 10 tanto de PC1 como de PC2
library(patchwork)
# En un grid 2x1
plot_pc1 + plot_pc2 + plot_layout(ncol = 2, nrow = 1)
GRÁFICAS PCA
Utilizamos el paquete factoextra para realizar las
gráficas de los 2 primeros componentes, en concreto cla función
fviz_pca_ind, que además nos permite añadir las elipses de
confianza.
Entre los diferentes grupos de estudio, las elipses de confianza están centradas y solapadas lo que sugiere que, en el espacio de los primeros dos componentes principales, no hay una separación clara entre los grupos.
No es descartable la necesidad de examinar componentes adicionales o utilizar otras técnicas de análisis como métodos de Clustering.
# Utilizamos "factoextra" para hacer la gráfica de la PCA
library(factoextra)
# Hacemos la gráfica de los dos primeros componentes
fviz_pca_ind(pca_result,
col.ind = df_dataClin$study, # Colores de los puntos basados en los estudios
label = "none",
addEllipses = TRUE, # Añadir elipses de confianza
ellipse.type = "confidence", # Tipo de elipse, puede ser 'confidence' o 't'
legend.title = "Study",
palette = "jco", # Paleta de colores, 'jco' es solo una opción
title = "PCA Estudios")
Cuando representamos la PCA por inestabilidad de microsatélites vemos
que las dos elipses, las de MSI y las de MSS, están claramente separadas
lo que sugiere que ambos grupos tienen perfiles de expresión genética
distintos en las dimensiones capturadas por los componentes principales
(CP).
Como ya sabemos, estas diferencias son biológicamente relevantes; los
componentes que muestran esta separación probablemente están capturando
las variaciones genéticas ya conocidas de antemano.
# Hacemos la gráfica de los dos primeros componentes en relación a los grupos
fviz_pca_ind(pca_result,
col.ind = df_dataClin$msi_imputed, # Colores de los puntos basados en MSI/MSS
label = "none",
addEllipses = TRUE, # Añadir elipses de confianza
ellipse.type = "confidence", # Tipo de elipse, puede ser 'confidence' o 't'
legend.title = "MSI/MSS",
palette = "jco", # Paleta de colores, 'jco' es solo una opción
title = "PCA MSI/MSS")
Entre los dos grupos de estadio, II y III, las elipses de confianza están más o menos centradas y algo solapadas lo que sugiere que, no hay una separación clara entre los grupos. Algo un tanto previsible ya que el estadio del cáncer depende del tiempo de progreso de la enfermedad.
# Hacemos la gráfica de los dos primeros componentes en relación a los grupos
fviz_pca_ind(pca_result,
col.ind = df_dataClin$stage, # Colores de los puntos basados en Estadio
label = "none",
addEllipses = TRUE, # Añadir elipses de confianza
ellipse.type = "confidence", # Tipo de elipse, puede ser 'confidence' o 't'
legend.title = "Estadio",
palette = "jco", # Paleta de colores, 'jco' es solo una opción
title = "PCA Estadio II/III")
Se podría decir que esta gráfica és una combinación de las dos anteriores; fijándonos en las elipses de confianza, se observan dos grupos separados claramente de otros dos, debido principalmente a la contribución de la instabilidad de microsatélites y no tanto a la del estadio del CRC.
# Hacemos la gráfica de los dos primeros componentes en relación a los grupos
fviz_pca_ind(pca_result,
col.ind = df_dataClin$grupo, # Colores de los puntos basados en Grupos
label = "none",
addEllipses = TRUE, # Añadir elipses de confianza
ellipse.type = "confidence", # Tipo de elipse, puede ser 'confidence' o 't'
legend.title = "Grupo",
palette = "jco", # Paleta de colores, 'jco' es solo una opción
title = "PCA Grupos")
Al igual que con la inestabilidad de microsatélites, cuando
representamos la PCA por clasificación CMS vemos que las cuatro elipses
están claramente separadas; los 4 grupos tienen perfiles de expresión
genética distintos en las dimensiones capturadas por los CP.
Estas diferencias son biológicamente relevantes dado que los componentes
que muestran esta separación están capturando variaciones genéticas ya
conocidas.
# Hacemos la gráfica de los dos primeros componentes en relación a los grupos
fviz_pca_ind(pca_result,
col.ind = df_dataClin$cms, # Colores de los puntos basados en CMS
label = "none",
addEllipses = TRUE, # Añadir elipses de confianza
ellipse.type = "confidence", # Tipo de elipse, puede ser 'confidence' o 't'
legend.title = "CMS",
palette = "jco", # Paleta de colores, 'jco' es solo una opción
title = "PCA por CMS")
Analizamos los datos viendo la infiltración celular.
Utilizamos el paquete MCPcounter (Becht et al. 2016) que estima la abundancia poblacional de poblaciones de células inmunes (8 tipos) y estromales (2 tipos) que se infiltran en tejidos mediante la expresión genética.
En concreto, utilizatemos MCPcounter para estimar la infiltración celular del estroma y el sistema inmunitario en el microambiente tumoral (TME).
# Cargamos la librería 'MCPcounter'
require(MCPcounter)
# Con todos los datos de la matriz de expresión 26/03/2024
exmat_MCP <- MCPcounter.estimate(expression = df_exmat, featuresType = "HUGO_symbols")
# Transponemos
exmat_MCP <- t_df(exmat_MCP)
# Añadimos los datos clínicos a los datos de la matriz generada con MCPcounter
# 26/03/2024
exmat_MCP_Clin <- cbind(df_dataClin, exmat_MCP)
A modo de recordatorio didáctico de Inferencia Estadística:
| \(H_0\) es verdadera | \(H_0\) es falsa | |
|---|---|---|
| Rechazar \(H_0\) | Error Tipo I (α) | Decisión Correcta |
| No rechazar \(H_0\) | Decisión Correcta | Error Tipo II (β) |
El p-valor es la probabilidad del resultado del estadístico de contraste observado o de uno más alejado cuando la hipótesis nula es cierta.
El p-valor asociado a una observación del estadístico de contraste es el menor nivel de significación que nos permite rechazar la hipótesis nula.
Si el p-valor es inferior al nivel de significación \(\alpha\), rechazaremos la hipótesis nula, en caso contrario la aceptaremos.
En nuestro caso utilizamos un nivel de significación \(\alpha = 0.05\)
Para hacer una comparativa entre los 4 grupos, utilizaremos el test no paramétrico de Kruskal-Wallis. Se suele utilizar para determinar si hay diferencias significativas entre tres o más grupos independientes cuando los datos no cumplen los supuestos necesarios para realizar un ANOVA: normalidad y homocedasticidad (homogeneidad de varianzas).
# Calculamos los test no paraméticos Kruskal-Wallis (asumimos que no tenemos
# normalidad) 16/04/2024
# tests de Kruskal-Wallis para NK
kw_NK_result <- kruskal.test(exmat_MCP_Clin$`NK cells` ~ grupo, data = exmat_MCP_Clin)
# Para el resto de tipos celulares utilizamos una lista con los nombres
tipos_celulares <- c("NK cells", "T cells", "CD8 T cells", "Cytotoxic lymphocytes",
"B lineage", "Monocytic lineage", "Myeloid dendritic cells", "Neutrophils", "Endothelial cells",
"Fibroblasts")
# Aplicamos el test de Kruskal-Wallis a cada tipo celular con 'lapply'
resultados_kw <- lapply(tipos_celulares, function(cell_type) {
kruskal.test(reformulate("grupo", response = cell_type), data = exmat_MCP_Clin)
})
# Añadimos los nombres de los tipos celulares
names(resultados_kw) <- tipos_celulares
# Extraemos los valores de p de todos los tests realizados
p_valores <- sapply(resultados_kw, function(x) x$p.value)
La hipótesis nula del test de Kruskal-Wallis \(H_0\) establece que todas las poblaciones (grupos) tienen distribuciones idénticas, es decir, las medianas de todos los grupos son iguales.
\[ H_0: \forall i, j \, (i \neq j) \Rightarrow \mu_i = \mu_j \]
La hipótesis alternativa \(H_1\) establece que al menos uno de los grupos tiene una distribución diferente en cuanto a la mediana, comparada con las otras poblaciones sin especificar cuál de los grupos es diferente.
\[ H_1: \exists i \neq j \mid \mu_i \neq \mu_j \]
# Generamos todos los boxplot 27/03/2024
grupo <- exmat_MCP_Clin$grupo
# NK
bp1 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`NK cells`, grupo,
"NK cells", p_valores["NK cells"])
# T cells
bp2 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`T cells`, grupo,
"T cells", p_valores["T cells"])
# CD8 T cells
bp3 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`CD8 T cells`, grupo,
"CD8 T cells", p_valores["CD8 T cells"])
# Cytotoxic lymphocytes
bp4 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`Cytotoxic lymphocytes`,
grupo, "Cytotoxic lymphocytes", p_valores["Cytotoxic lymphocytes"])
# B lineage
bp5 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`B lineage`, grupo,
"B lineage", p_valores["B lineage"])
# Monocytic lineage
bp6 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`Monocytic lineage`,
grupo, "Monocytic lineage", p_valores["Monocytic lineage"])
# Myeloid dendritic cells
bp7 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`Myeloid dendritic cells`,
grupo, "Myeloid dendritic cells", p_valores["Myeloid dendritic cells"])
# Neutrophils
bp8 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$Neutrophils, grupo,
"Neutrophils", p_valores["Neutrophils"])
# Endothelial cells
bp9 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$`Endothelial cells`,
grupo, "Endothelial cells", p_valores["Endothelial cells"])
# Fibroblasts
bp10 <- boxplot_MCPcounter(exmat_MCP_Clin, grupo, exmat_MCP_Clin$Fibroblasts, grupo,
"Fibroblasts", p_valores["Fibroblasts"])
library(patchwork)
# En un grid 3x2
bp1 + bp2 + bp3 + plot_layout(ncol = 3, nrow = 1)
bp4 + bp5 + bp6 + plot_layout(ncol = 3, nrow = 1)
bp7 + bp8 + bp9 + plot_layout(ncol = 3, nrow = 1)
bp10 + plot_layout(ncol = 3, nrow = 1)
NORMALIDAD Y HOMOCEDASTICIDAD
Primero compromabos la normalidad de los grupos para las NK cells.
Tabla tabla con los resultados de los tests de
Shapiro para normalidad.
Si p-value < 0.05 rechazamos la hipótesis nula de igualdad de
estadísticos.
| Shapiro test | MSI | MSS | II | III |
|---|---|---|---|---|
| p-valor | 3.339e-05 | 1.058e-10 | 1.371e-10 | 1.089e-04 |
Dado que no hay normalidad no hace falta comprobar homocedasticidad.
Para realizar una comparativa entre dos grupos en los que asumimos no normalidad , en nuestro caso MSI vs MSS y estadio II vs estadio III, utilizamos el test no paramétrico de Mann-Whitney U también conocido como el test de Wilcoxon en el que se asume que los datos no necesitan ser normalmente distribuidos aunque deben poder ordenarse.
# Para dos grupos utilizamos el test no paramético Mann-Whitney U (asumimos que
# no tenemos normalidad) 16/04/2024
# Realizamos el test de Mann-Whitney U para los grupos MSI y MSS
mwU_NK_msi_test <- wilcox.test(grupo_msi, grupo_mss, alternative = "two.sided")
# Realizamos el test de Mann-Whitney U para los estadios II y III
mwU_NK_stage_test <- wilcox.test(grupo_sII, grupo_sIII, alternative = "two.sided")
La hipótesis nula \(H_0\) para el test de Mann-Whitney U establece que no hay diferencia en las medianas (o, más generalmente, en las distribuciones) de los dos grupos: \[ H_0: \mu_{grupo1} = \mu_{grupo2} \] La hipótesis alternativa \(H_1\) establece que hay una diferencia entre las distribuciones de los dos grupos. En nuestro caso la planteamos Bilateral (Two-sided):
\[ H_1: \mu_{grupo1} \neq \mu_{grupo2} \]
# Generamos todos los boxplot referentes a NKs 27/03/2024 Diagrama de Barras de
# NKs MSI vs MSS
bp11 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`NK cells`,
exmat_MCP_Clin$msi_imputed, "NK cells MSI vs MSS", mwU_NK_msi_test$p.value)
# Diagrama de Barras de NKs II vs III
bp12 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`NK cells`,
exmat_MCP_Clin$stage, "NK cells II vs III", mwU_NK_stage_test$p.value)
require(patchwork)
# En un grid 2x2
p4 + bp1 + p1 + bp11 + p3 + bp12 + plot_layout(ncol = 2, nrow = 3)
# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp13 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`T cells`,
exmat_MCP_Clin$msi_imputed, "T cells MSI vs MSS")
# Diagrama de Barras II vs III
bp14 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`T cells`,
exmat_MCP_Clin$stage, "T cells II vs III")
require(patchwork)
# En un grid 1x3
bp2 + bp13 + bp14 + plot_layout(ncol = 3, nrow = 1)
# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp15 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`CD8 T cells`,
exmat_MCP_Clin$msi_imputed, "CD8 T cells MSI vs MSS")
# Diagrama de Barras II vs III
bp16 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`CD8 T cells`,
exmat_MCP_Clin$stage, "CD8 T cells II vs III")
require(patchwork)
# En un grid 1x3
bp3 + bp15 + bp16 + plot_layout(ncol = 3, nrow = 1)
# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp17 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`Cytotoxic lymphocytes`,
exmat_MCP_Clin$msi_imputed, "Cytotoxic lymphocytes MSI vs MSS")
# Diagrama de Barras II vs III
bp18 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`Cytotoxic lymphocytes`,
exmat_MCP_Clin$stage, "Cytotoxic lymphocytes II vs III")
require(patchwork)
# En un grid 1x3
bp4 + bp17 + bp18 + plot_layout(ncol = 3, nrow = 1)
# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp19 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`B lineage`,
exmat_MCP_Clin$msi_imputed, "B lineage MSI vs MSS")
# Diagrama de Barras II vs III
bp20 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`B lineage`,
exmat_MCP_Clin$stage, "B lineage II vs III")
require(patchwork)
# En un grid 1x3
bp5 + bp19 + bp20 + plot_layout(ncol = 3, nrow = 1)
# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp21 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`Monocytic lineage`,
exmat_MCP_Clin$msi_imputed, "Monocytic lineage MSI vs MSS")
# Diagrama de Barras II vs III
bp22 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`Monocytic lineage`,
exmat_MCP_Clin$stage, "Monocytic lineage II vs III")
require(patchwork)
# En un grid 1x3
bp6 + bp21 + bp22 + plot_layout(ncol = 3, nrow = 1)
# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp23 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`Myeloid dendritic cells`,
exmat_MCP_Clin$msi_imputed, "Myeloid dendritic cells MSI vs MSS")
# Diagrama de Barras II vs III
bp24 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`Myeloid dendritic cells`,
exmat_MCP_Clin$stage, "Myeloid dendritic cells II vs III")
require(patchwork)
# En un grid 1x3
bp7 + bp23 + bp24 + plot_layout(ncol = 3, nrow = 1)
# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp25 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$Neutrophils,
exmat_MCP_Clin$msi_imputed, "Neutrophils MSI vs MSS")
# Diagrama de Barras II vs III
bp26 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$Neutrophils,
exmat_MCP_Clin$stage, "Neutrophils II vs III")
require(patchwork)
# En un grid 1x3
bp8 + bp25 + bp26 + plot_layout(ncol = 3, nrow = 1)
# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp27 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$`Endothelial cells`,
exmat_MCP_Clin$msi_imputed, "Endothelial cells MSI vs MSS")
# Diagrama de Barras II vs III
bp28 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$`Endothelial cells`,
exmat_MCP_Clin$stage, "Endothelial cells II vs III")
require(patchwork)
# En un grid 1x3
bp9 + bp27 + bp28 + plot_layout(ncol = 3, nrow = 1)
# Generamos todos los boxplot referentes 01/04/2024 Diagrama de Barras MSI vs
# MSS
bp29 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$msi_imputed, exmat_MCP_Clin$Fibroblasts,
exmat_MCP_Clin$msi_imputed, "Fibroblasts MSI vs MSS")
# Diagrama de Barras II vs III
bp30 <- boxplot_MCPcounter(exmat_MCP_Clin, exmat_MCP_Clin$stage, exmat_MCP_Clin$Fibroblasts,
exmat_MCP_Clin$stage, "Fibroblasts II vs III")
require(patchwork)
# En un grid 1x3
bp10 + bp29 + bp30 + plot_layout(ncol = 3, nrow = 1)
# En un grid 1x3
bp1 + bp11 + bp12 + plot_layout(ncol = 3, nrow = 1)
La clasificación CMS divide los tumores de CRC en cuatro subtipos principales, basados en características moleculares distintas, lo que refleja diferencias en la biología del tumor, la interacción con el TME, el pronóstico y la respuesta a tratamientos específicos. (Roelands et al. 2017)
| Subtipo | Características | Inmunidad | Prevalencia |
|---|---|---|---|
| CMS1 (Subtipo Inmuno) | Alta inmunogenicidad, frecuentemente asociada con MSI y una alta carga mutacional. | Tienen una infiltración activa de células inmunitarias efectoras, lo que puede hacerlos más susceptibles a terapias inmunogénicas. | Aproximadamente el 14% de los casos de CRC. |
| CMS2 (Subtipo Canónico) | Tumores con alta expresión de genes relacionados con el ciclo celular y la señalización WNT/MYC, presentando una marcada proliferación celular. | Son generalmente poco inmunogénicos, con menor infiltración inmune comparado con CMS1 y CMS4. | Es el subtipo más común, abarcando alrededor del 37% de los casos de CRC. |
| CMS3 (Subtipo Metabólico) | Presenta alteraciones metabólicas, incluyendo un metabolismo energético distintivo, con una menor carga mutacional en comparación con CMS1. | Similar a CMS2, muestra baja inmunogenicidad. | Constituye aproximadamente el 13% de los casos de CRC |
| CMS4 (Subtipo Mesenquimal) | Se caracteriza por la activación de vías relacionadas con la angiogénesis, la activación de células estromales y la señalización TGF-β, lo que conduce a un entorno inmunosupresor. | Aunque presenta una significativa infiltración inmune, esta es predominantemente de naturaleza supresora, con una alta presencia de células estromales. | Aproximadamente el 23% de los casos de CRC caen en este subtipo. |
Vemos la infiltración celular según CMS.
# Calculamos los test no paraméticos Kruskal-Wallis (asumimos que no tenemos
# normalidad) 16/04/2024
# Aplicamos el test de Kruskal-Wallis a cada tipo celular con 'lapply'
resultados_kw_cms <- lapply(tipos_celulares, function(cell_type) {
kruskal.test(reformulate("cms", response = cell_type), data = exmat_MCP_Clin)
})
# Añadimos los nombres de los tipos celulares
names(resultados_kw_cms) <- tipos_celulares
# Extraemos los valores de p de todos los tests realizados
p_valores_cms <- sapply(resultados_kw_cms, function(x) x$p.value)
# Quitamos los NA's de la columna cms
exmat_MCP_Clin_CMS <- exmat_MCP_Clin[!is.na(exmat_MCP_Clin$cms), ]
# Generamos todos los boxplot 02/04/2024
cms <- exmat_MCP_Clin_CMS$cms
# NK
bp1 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`NK cells`,
cms, "NK cells", p_valores_cms["NK cells"])
# T cells
bp2 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`T cells`,
cms, "T cells", p_valores_cms["T cells"])
# CD8 T cells
bp3 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`CD8 T cells`,
cms, "CD8 T cells", p_valores_cms["CD8 T cells"])
# Cytotoxic lymphocytes
bp4 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`Cytotoxic lymphocytes`,
cms, "Cytotoxic lymphocytes", p_valores_cms["Cytotoxic lymphocytes"])
# B lineage
bp5 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`B lineage`,
cms, "B lineage", p_valores_cms["B lineage"])
# Monocytic lineage
bp6 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`Monocytic lineage`,
cms, "Monocytic lineage", p_valores_cms["Monocytic lineage"])
# Myeloid dendritic cells
bp7 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`Myeloid dendritic cells`,
cms, "Myeloid dendritic cells", p_valores_cms["Myeloid dendritic cells"])
# Neutrophils
bp8 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$Neutrophils,
cms, "Neutrophils", p_valores_cms["Neutrophils"])
# Endothelial cells
bp9 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$`Endothelial cells`,
cms, "Endothelial cells", p_valores_cms["Endothelial cells"])
# Fibroblasts
bp10 <- boxplot_MCPcounter(exmat_MCP_Clin_CMS, cms, exmat_MCP_Clin_CMS$Fibroblasts,
cms, "Fibroblasts", p_valores_cms["Fibroblasts"])
require(patchwork)
# En un grid 3x2
p2 + bp1 + bp2 + plot_layout(ncol = 3, nrow = 1)
bp3 + bp4 + bp5 + plot_layout(ncol = 3, nrow = 1)
bp6 + bp7 + bp8 + plot_layout(ncol = 3, nrow = 1)
bp9 + bp10 + plot_layout(ncol = 3, nrow = 1)
Los tumores CMS2 y CMS3 son poco inmunogénicos; CMS1 y CMS4 difieren
en el tipo de infiltración inmune.
Los tumores CMS1 tienden a acumular una gran cantidad de mutaciones
debido al estado de MSI y atraen células efectoras inmunes.
Los tumores CMS4 exhiben un TME inmunosupresor con infiltración
estromal. (Lanuza et al. 2022)
Repositorio: (Ferrón 2024)